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Abstract. We review the properties of reduced density matrices for free fermionic or 
bosonic many-particle systems in their ground state. Their basic feature is that they 
have a thermal form and thus lead to a quasi-thermodynamic problem with a certain 
C/2 . free-particle Hamiltonian. We discuss the derivation of this result, the character of the 

Hamiltonian and its eigenstates, the single-particle spectra and the full spectra, the 
resulting entanglement and in particular the entanglement entropy. This is done for 
various one- and two-dimensional situations, including also the evolution after global 
or local quenches. 

o 
o 

^ ! 1. Introduction 

> . 

Reduced density matrices contain the information on some part of a quantum system 
and are a basic tool in many-body physics. The ones commonly employed describe 
the properties of one or two selected particles in a many-particle system and allow to 
O ' calculate important physical quantities like the total energy or the density correlations. 
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These reduced density matrices (RDM's) were first introduced by Dirac pQ and studied 
already in the 1930's, see e.g. [2, E]. In the usual terminology, they are just the static 
one- and two-particle correlation functions. 

The RDM's we want to discuss here are of a different type and refer to a different 
question. They arise if one divides a system in space, or, more generally, in Hilbert space, 
and asks how the two parts are coupled in the given wave function. This corresponds to 
the analysis by Schrodinger in 1935 [I] when he introduced the concept of entanglement. 
The general form of this coupling is given by the Schmidt decomposition which displays 
all entanglement features in a simple and transparent way. To obtain it in a specific 
case, one needs the RDM's for the two regions in question. 

The present interest in this problem, although it had also been a theme in quantum 
optics, arose in the beginning of the nineties in two seemingly disconnected areas, in the 
theory of black holes [3 [6], [7] and in the numerical investigation of quantum chains [H[9]. 
In both cases, the motivation came from the wish to consider some subsystem which 
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is in contact with its environment. For the quantum chains, this lead to the density- 
matrix renormalization group (DMRG) which can treat large systems with spectacular 
accuracy and revolutionized the field [TQ(, [TT]. A third input then came from the area 
of quantum information, where the structure of quantum states also plays a central 
role. This resulted in particular in a renewed and extensive study of the entanglement 
entropy [T2"j [T3| [H] which is a simple and convenient measure of the entanglement and 
follows directly from the RDM eigenvalues. 

The purpose of this article is to give a coherent account of the reduced density 
matrices just described for a class of models where they can be obtained in closed 
form. These are free fermions including the related spin chains and free bosons in the 
form of coupled oscillators. They will be considered either in their ground state or 
in certain other pure quantum states. In this case, the RDM's are found to have a 
Boltzmann-like form with a certain free-particle operator in the exponent. The problem 
is thereby reduced to the study of this associated Hamiltonian and its characteristic 
features. The main property of interest is the eigenvalue spectrum since it determines 
the spectrum of the RDM itself and thus the entanglement properties, in particular the 
entanglement entropy. Both the spectra and the entropies will be presented for a variety 
of different situations. The problem on a lattice is very clear-cut. The partitioning is 
done by selecting two sets of discrete sites and there are no divergencies for finite sizes. 
On the other hand, there is only a small number of analytical results and one has to 
invoke numerics frequently. Lattice systems are also required for the DMRG, and the 
initial motivation for the studies was to understand the performance of this intriguing 
numerical method by looking at solvable models. 

In section 2 we will provide some background on entanglement, the Schmidt 
decomposition and the RDM's. Then, in section 3, we give the general form of the 
reduced density matrices for free fermions or bosons and discuss the methods for 
obtaining them. For quantum chains, this also contains relations to two-dimensional 
classical models. In section 4 we show the eigenvalue spectra for various one- and two- 
dimensional systems and discuss their typical appearance, their scaling behaviour and 
the change with the dimension. The characteristics of the single-particle eigenfunctions, 
the nature of the effective Hamiltonian and some further aspects are the topics of section 
5. In section 6 we turn to the entanglement entropy and summarize the important 
results with emphasis on their relation to the spectra. Finally, in section 7, we review 
the temporal behaviour of the entanglement after different types of quenches. The 
material is drawn preferentially from our own studies and some of it also appeared in 
a recent book [15]. However, the scope is different here and a considerable number of 
figures were prepared exclusively for this review. 

2. Background 

In this section, we summarize the basic features of entangled states and reduced density 
matrices in order to create the frame for the results to be presented later. For more 
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details, see e.g. the short review [IE] . 



2.1. Schmidt decomposition 

Consider a quantum system which is divided into two distinct parts 1 and 2. Then a 
state of the total system can be written 

m,n 

where l^ 1 .) an d \^n) are orthonormal basis functions in the two Hilbert spaces. But 
a rectangular matrix A can always be written in the form UDV where U is unitary, 
D is diagonal and the rows of V are orthonormal. This is called the singular- value 
decomposition and similar to the principal-axis transformation of a symmetric square 
matrix [17J. Using this in ([1]) and forming new bases by combining the |^/^) with U 
and the with V, one obtains the Schmidt decomposition [TB] 

I^HE^I^)!^) ( 2 ) 

n 

which gives the total wave function as a single sum of products of orthonormal functions. 
Here the number of terms is limited by the smaller of the two Hilbert spaces and the 
weight factors A n are the elements of the diagonal matrix D. If is normalized, their 
absolute magnitudes squared sum to one. The entanglement properties are encoded in 
the set of A n . Only if all except one are zero, the sum reduces to a single term and |^) is 
a product state, i.e. non-entangled. In all other cases a certain entanglement is present 
and if all \ n are equal in size, one would call the state maximally entangled. Of course, 
this refers to a particular bipartition and one should investigate different partitions to 
obtain a complete picture. 



2.2. Reduced density matrices 

The entanglement structure just discussed can also be found from the density matrices 
associated with the state This is, in fact, the standard way to obtain it. Starting 
from the total density matrix 

P=I*X*I (3) 

one can, for a chosen division, take the trace over the degrees of freedom in one part of 
the system. This gives the reduced density matrix for the other part, i.e. 

Pi = tr 2 (p) , p 2 = tix(p) (4) 

These hermitian operators can be used to calculate arbitrary expectation values in the 
subsystems. Moreover, it follows from (T5]) that their diagonal forms are 

Pa = Y,\ X n\ 2 \®n)(K\ , OC = 1,2 (5) 

n 

This means that 



Reduced density matrices and entanglement entropy in free lattice models 



4 



• pi and P2 have the same non-zero eigenvalues 

• these eigenvalues are given by w n = |A n | 2 

Therefore the eigenvalue spectrum of the p a gives directly the weights in the Schmidt 
decomposition and a glance at this spectrum shows the basic entanglement features 
of the state, for the chosen bipartition. For this reason, it has also been termed 
"entanglement spectrum" recently [19] . One also sees that the |$") appearing in 
([2]) are the eigenfunctions of the p a . For the single-particle RDM's mentioned in the 
introduction, these eigenfunctions are known as "natural orbitals" in quantum chemistry 



In the DMRG algorithm, these properties are used to truncate the Hilbert space 
by calculating the p a , selecting the m states |$") with largest weights w n and deleting 
the rest. This procedure is expected to work well if the total weight of the discarded 
states is sufficiently small. Therefore the form of the density-matrix spectra is decisive 
for the success of the method. 

It is interesting that Schmidt himself already worked with the RDM's. Studying 
coupled linear integral equations, he derived a spectral representation of the form (T5]) for 
an unsymmetric kernel K in terms of the eigenfunctions of the two symmetric operators 
KK' and K'K. His paper (which is based on his doctoral thesis with Hilbert) also 
contains the recipe for the best approximation as it is used in the DMRG. 

2.3. Entanglement entropy 

Whereas the full RDM spectra give the clearest impression of the entanglement in a 
bipartite system, it is also desirable to have a simple measure which condenses this 
information into one number. This can be achieved by generalizing the usual (von 
Neumann) entropy definition to reduced density matrices. The entanglement entropy 
therefore reads: 



where the trace has been rewritten as a sum using the eigenvalues w n . The most 
important properties are as follows. 

• The entropy is determined purely by the spectrum of pi, which is known to be 
identical to the spectrum of p^. Therefore S\ = 5*2 holds for arbitrary bipartitions 
and one can simply write S and talk of the entanglement entropy. 

• The entropy vanishes for product states, and has a maximal value of S = In M if 
one has M non-zero eigenvalues which are all equal, w n = 1/M for n — 1,2, ... , M. 
Using this, one can write in general S = lnM e g, thereby defining an effective 
number of states coupled in parts 1 and 2. This gives a simple interpretation to S. 

Although there are other entanglement measures [21] , the entropy is the standard 
one for bipartitions and will be discussed in detail later. It is important to keep in mind 
that it measures a mutual connection and will, in general, not be proportional to the 
size of a subsystem. 



1201. 




(6) 



n 
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3. RDM's for free lattice models 



3.1. Systems 



In the following we consider models with a Hamiltonian which is quadratic in 
either fermion or boson operators and thus can be diagonalized by a Bogoliubov 
transformation. In principle, these can be quite general, but we will concentrate on 
the following physically important systems 

• Fermionic hopping models with conserved particle number and Hamiltonian 

H = ~ 2 ^ ] t m ^ n (? m c n (7) 

<m,n> 

where the symbol <> denotes nearest neighbours. Apart from homogeneous 
systems we will consider dimerized chains, where £ n ,n+i alternates between 1 ± 8, 
and the case of single defects. 

• Coupled oscillators with eigenfrequency ujq and Hamiltonian 



H 



^ ~9 + 7j u o x l + T k m ,n(x m -x n ) 2 (8) 



2dxl 2 



4 

<m,n> 



These are systems with an optical spectrum and bosonic pair creation and 
annihilation. 

Spin one-half chains which are equivalent to free fermions via the Jordan- Wigner 
transformation. The most general one is the XY chain with a Z field, described by 



7 „x„x i T _« „y 

-T—Vn a n+1 + — <°n+l 



hJ2< ( 9 ) 



where the cr" are Pauli matrices at site n. For 7 = this reduces to the XX model, 
corresponds to ((?]) with nearest-neighbour hopping and can also model hard-core 
bosons. For 7 7^ 0, it contains pair creation and annihilation terms. For 7 = 1 
it becomes the Ising model in a transverse field (TI model) which we write, in a 
slightly different notation 

n n 

The solubility of the models in itself does not yet mean that the RDM's are easily 
accessible. For example, they have been considered in the critical XXZ spin chain, but 
the formulae are very complicated, see [221 123] . The free lattice models, however, have 
eigenstates with special properties which permit to make a simple general statement. 

3.2. General result 

For these free-particle models, the reduced density matrices for the ground state can be 
written 

1 L 
Pa = ^e- n -, H a = J2 £ if!fi (U) 

1=1 
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Here L is the number of sites in subsystem a and the operators fj ,fi are fermionic or 
bosonic creation and annihilation operators for single-particle states with eigenvalues 
E\. The /'s are related to the original operators in the subsystem by a canonical 
transformation. Thus p a has the form of a thermal density matrix with an effective 
Hamiltonian 7i a which is of the same free-particle type as H. In ( TTTT) it is already given 
in diagonal form. The constant Z, written in analogy to thermodynamics, ensures the 
correct normalization tr(p a ) = 1. 

This form of p a is rather suggestive since one has a similar situation as for a 
system in contact with a thermal bath. However, no assumption about the relative 
sizes of the two coupled systems enters here. More importantly, the operator 7i a is 
not the Hamiltonian H restricted to the subsystem a. Therefore (TIT]) is not a true 
Boltzmann formula. Nevertheless, the problem has been reduced to the study of a 
certain Hamiltonian and its thermodynamic properties. The features of l~L a will be 
the topic of the next chapters. Generally, one can say that it corresponds to an 
inhomogeneous system even if the subsystem it describes is homogeneous. This will 
be seen in more detail in section 5.2. Here we first discuss how one arrives at ffTTj) . 
These considerations will also show that validity of ( TTTT) goes even beyond the ground 
state. 

3.3. Methods 

Basically, there are three methods to obtain the reduced density matrices. 

(I) Integration over part of the variables according to the definition (jl]). This can 
be done e.g. for N coupled harmonic oscillators [2^, [25]. In this case the ground state is 
a Gaussian in the normal coordinates, provided no normal frequency vanishes. In terms 
of the original coordinates x n of the oscillators, it has the form 



Here C is a normalization constant and the matrix A is the square root V 1 / 2 of the 
dynamical matrix associated with the potential energy. By forming p and integrating 
out e.g. the variables xl+i, . . . ,xn one obtains p%(xi, X2, ■ ■ ■ , xl | x[, x' 2 , ■ ■ ■ , x' L ) which 
is again a Gaussian. With proper linear combinations yi of the coordinates, it contains 
only squares yf, yf and differences (yi — y[) 2 . Early treatments worked with this integral 
operator [5j E]. However, one can convert the differences into second derivatives and 
thereby obtain the differential operator 



where the exponents become quadratic expressions in terms of boson operators. A 
diagonalization then gives the single exponential (1111) with Tii describing a collection of 
L new harmonic oscillators. Their eigenfrequencies Si follow from A by dividing it into 




(12) 




L 



(13) 
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the submatrices a 11 , a 12 , a 21 , and a 22 , according to whether the sites are in part 1 or in 
part 2. Then the LxL matrix a 11 (a 11 — a 12 (a 22 )^ 1 a 21 )^ 1 has the eigenvalues coth 2 (e//2). 

If L = 1, there is just one such oscillator with a frequency e which differs from ujq. 
Its eigenstates have a different spatial extent and may therefore be called "squeezed". 
For N = 2 the resulting Schmidt decomposition of \I> (xi, x%) in terms of these states can 
easily be written down and is well known, see e.g. [26j [27] . 

The method can also be used for systems of non-interacting fermions. In this case 
one first has to write the ground state in exponential form and then use Grassmann 
variables for the integration [28l [29] . 

(II) Via correlation functions [30~1 f!2l l3T] . The simplest case is a system of free 
electrons hopping on TV lattice sites in a state described by a Slater determinant. In such 
a state, all many-particle correlation functions factorize into products of one-particle 
functions. For example, 

(c^cickCi) = {c ] m ci){c ] n c k ) - (c ] m c k )(c ] n d) (14) 

If all sites are in the same subsystem, a calculation using the reduced density matrix 
must give the same result. This is guaranteed by Wick's theorem if p a is the exponential 
of a free-fermion operator 

L 

p a = Kexp (- ^2 Ki c \ c i) ( 15 ) 

where % and j are sites in the subsystem. With the form of p a fixed, the hopping matrix 
hij is then determined such that it gives the correct one-particle correlation functions 
Cij = (cjcj). The two matrices are diagonalized by the same transformation and one 
finds (see also [29 J) 

h = ln [(1-C)/C] (16) 

The same formula also relates the eigenvalues E\ and (i of h and C. Expressed differently, 
the Ei follow from the equation 

(l-2C)0j = tanh(|)0*. (17) 

If there is pair creation and annihilation, one has to include the 'anomalous' 
correlation functions Fij = (cjcj-) and F*j = (cjCi). To reproduce them, the operator 
7i a then must also contain pair terms. Diagonalizing it in the usual way [32], one finds 
that the single-particle eigenvalues follow from two coupled equations, which can be 
combined into a single one. For real F this reads 

(2C - 1 - 2F)(2C - 1 + 2F) ; = tanh 2 (^) 0,. (18) 

and reduces to the previous result (fTTl) if F vanishes. Alternatively, one can work 
with Majorana operators [121 [13] &2n-i = (c n + cjj and a<m = i(c n — cjj and form 
the 2N x 2N correlation matrix M m>n = (a m a n ). Restricted to the subsystem, it 
contains the same elements as the two matrices in (jT8j) but arranged differently. Writing 
M mtn = 5 myn + iT m ^ n , the matrix T of the subsystem has the eigenvalues ±i tanh(^/2). 
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This method is very general. It works in any dimension, for arbitrary quadratic 
Hamiltonians, for all states which are Slater determinants, and even at finite 
temperature. Thus it has been used in a large number of situations ranging from 
homogeneous chains to defect problems, random systems, higher dimensions and the 
time evolution after a quench. 

Factorization properties as in (j!4p are well-known for Gaussians, and therefore the 
approach is equally applicable to coupled oscillators in the ground state f|T2|) . Thus 
p a must be the exponential of a bosonic operator (as found in (I)) and 7i a is again 
determined such that it reproduces the correlation functions, in this case those of 
positions and momenta, X^j = (xiXj) and Pij = (piPj). In analogy to ( fTSl the single- 
particle eigenvalues then follow from [33], [301 EI] 

2P2X^ = coth 2 (|)0 ; . (19) 

Since for the total system 2V = V 1 / 2 = A and 2X = V -1 / 2 , the matrix on the left side 
of (|T9|) is seen to be the restriction of A to the subsystem multiplied by the restriction of 
its inverse. This is exactly the expression given in (I). As in the fermionic case, one can 
also combine coordinates and momenta, which are analogous to the Majorana variables, 
and consider the corresponding 2L x 2L correlation matrix, usually called covariance 
matrix. Its reduction to diagonal form is a well-known problem in mathematics [35j and 
the resulting coth(e//2) are also referred to as symplectic eigenvalues [36l 137]. 

The method was used for example in [3B1 [39], [31] and again works also at finite 
temperature. 

(Ill) Via classical statistical models [40] [41]. In one dimension one can exploit the 
relations between quantum chains and two-dimensional classical models. The starting 
point is a discrete version of a path-integral representation. 

Consider a quantum chain of finite length and imagine that one can obtain its 
state |^) from an initial state \ty s ) by applying a proper operator T many times. If 
T is the row-to-row transfer matrix of a classical model, one has thereby related 
to the partition function of a two-dimensional semi-infinite strip of that system. The 
total density matrix is then given by two such strips. This is sketched on the 

far left of Fig{T] The reduced density matrix, e.g. for the left part of the chain, follows 
by identifying the variables along the right part of the horizontal edges and summing 
them, which means tying the two half-strips together. In this way, p a is expressed as 
the partition function of a full strip with a perpendicular cut, as shown half left in the 
figure. 

This procedure works for the ground state of a number of integrable quantum 
chains. For example, the TI chain can in this way be related to a two-dimensional Ising 
model on a square lattice which is rotated by 45° with respect to the horizontal |41j. 
In the same way, a chain of coupled oscillators is connected with a two-dimensional 
Gaussian model [21] and an XY chain with an Ising model on a triangular lattice [42J. 
Analogous correspondences link XXZ, XYZ and higher-spin chains to vertex models 
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Figure 1. Left: Density matrices for a quantum chain as two-dimensional partition 
functions. Far left: Expression for p. Half left: Expression for p\. The matrices 
are defined by the variables along the thick lines. Right: Two-dimensional system 
built from four quadrants with corresponding corner transfer matrices A, B, C, D. The 
arrows indicate the direction of transfer. After Ref. [15] . 



[4TI 1431 144] . To use these relations, however, one needs a way to actually calculate the 
resulting partition function. This is possible with the help of the corner transfer matrices 
(CTM's) introduced by Baxter [15]. These are partition functions of whole quadrants 
as shown on the right of FiglH or of sextants, if one is dealing with a triangular lattice. 
By multiplying these transfer matrices one can then obtain the reduced density matrix 
for a half-chain as 

p a ~ ABCD. (20) 

Since p a is given by an infinite strip, one also needs infinite-size CTM's in this relation. 
But exactly in this limit they are known for several non-critical integrable models and 
have the form 

A = e - uHc ™ (21) 

where u contains the anisotropy of the two-dimensional system. This is a consequence of 
the star-triangle relations on which the integrability rests [46] . This approach gives 7i a 
in the original variables, see section 5.2, and explicit expressions for the single-particle 
eigenvalues Si in the diagonalized form. According to the derivation, it applies to one- 
half of an infinite chain, but in practice the chain has only to be much longer than the 
correlation length. 

Summing up, we have shown how to arrive at (jTTI) and how to obtain the Si. The 
eigenstates of p a and their eigenvalues w n then follow by specifying the occupation 
numbers of all single-particle levels. The analytical result for Ei just mentioned is 
exceptional. For finite subsystems beyond one or two sites, one has to find the E\ 
numerically. This leads to a characteristic difficulty, because the eigenvalue equations in 
(II) contain hyperbolic functions which approach ±1 for large E\. As the subsystem size 
grows, more and more values lie (exponentially) close to ±1, and can only be obtained 
reliably with special techniques [47] . Therefore the values of the E\ in most of the 
following figures do not exceed 20-30. 
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4. Spectra 

In this section we give an overview of the single-particle spectra and the full p Q -spectra 
for various situations. These include different dimensions, critical and non-critical 
systems and the geometrical shape of the subsystem. We will focus on the Si because 
these are the primary quantities. 



4-1. One dimension 

(I) Non-critical chains. 

For infinite TI, XY and oscillator half-chains, the CTM approach gives the universal 
formula 



ei 



(2Z + l)e 
2le 



disordered region 
ordered region 



(22) 



where / = 0, 1, 2, . . .. Thus one has equidistant levels and in a plot E\ vs. I the dispersion 
is strictly linear. The only free parameter is the level spacing which depends on the 




Figure 2. Level spacing as a function of the parameter k. 



details of the model. It is given by 

e = 7rl{k')/l{k), (23) 

where I(k) denotes the complete elliptic integral of the first kind, and k' = \A — k 2 . 
The elliptic modulus k with < k < 1 is given in the TI model by 

In the XY model, the ordered region is subdivided by the so-called disorder line 
7 2 + h 2 = 1 and one has to distinguish three cases 



7 /y / 7 2 + fr 2 - l , h>l 
k = { v/ 7 2 + h 2 - I/7 , 7 2 + h 2 > 1, h < 1 (25) 

V(l - 7 2 -h 2 )/{\ -h 2 ) , 7 2 + /i 2 < l,h< 1 
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Here the last formula comes from a different approach [1BJ US]- For the oscillator chain, 
k is the nearest-neighbour coupling and one has to put ujq = 1 — k. In this case, there 
is no ordered region. In all models, the critical point is given by k = 1 and since I(k) 
diverges for k — > 1, the level spacing vanishes there and the dispersion curve becomes 
flat. The complete behaviour of e is shown in Fig. [2j 

Results for finite TI chains are shown in Fig. [3] on the left. The linear behaviour is 
perfect for the smallest A. As one comes closer to the critical point, the slope decreases 
as predicted, but there are also deviations from the linearity for large E\. Thus the linear 
region shrinks and is no longer visible at the critical point. This is the typical finite-size 
scenario in these models. On the right side, the resulting w n , ordered by magnitude, 
are shown. One can see a rapid decrease with n which is fastest for the smallest A but 
is still impressive at criticality (note the vertical scale). This means that a Schmidt 
decomposition could be truncated safely after about 10 terms and is the basis for the 
fantastic performance of the DMRG in this case [50J. 




1 2 3 4 5 6 7 8 9 10 10 20 30 40 



1 n 

Figure 3. Density-matrix spectra for one-half of a transverse Ising chain with N = 20 
sites in its ground state. Left: All ten single-particle eigenvalues £/. Right: The largest 
total eigenvalues w n - Reprinted with permission from [28]. ©2001 by the APS. 

The lowest w n -curve also shows a step structure with plateaus which become longer 
with n. These are a consequence of the equidistant levels, a certain eigenvalue of TC a can 
then be realized by different combinations of E\. The degeneracy is given by the number 
of partitions P(s) of an integer s into other (odd or even) integers. Using asymptotic 
formulae for the P{s), one finds the leading large-n behaviour [51] 

w n ~ exp[— a(lnn) 2 ] (26) 

where a = e6/tt 2 . The same result with a different constant a holds for bosons. If the 
dispersion is not strictly linear, the steps are smeared and a rather smooth w n spectrum 
is obtained. 

An important new feature appears in the ^-spectra, if the subsystem is a segment 
in a chain. Then a two-fold degeneracy is found, at least for the lowest eigenvalues. The 



Reduced density matrices and entanglement entropy in free lattice models 



12 



reason lies in the form of the eigenfunctions, which are concentrated near the ends, as 
will be demonstrated in section 5. This leads to a degeneracy of the w n , with a factor 
of 2 for each E\ which is involved, and therefore to a significantly slower decay. 

For the spin chains, there are cases where the ground state simplifies and becomes a 
doublet of product states. Then one Si is zero, while all others diverge. As a consequence, 
all w n except two collapse to zero. This happens not only in the TI model for A — > oo, 
but also in the XY model on the disorder line [28]. If the result were not known, one 
could locate the line from the behaviour of the spectra. 

Finally, we note that also a dimerized half-filled hopping model shows such 
equidistant e\ because one can relate it to the TI model via the correlation functions. 
The parameter k is then given by k = (1 — S)/(l + 5), where 5 > is the dimerization 
parameter. 

(II) Critical chains. 

In critical systems, the size of the subsystem affects not only the upper part of the 
single-particle spectrum. This is shown in Fig. H] for a segment in a half-filled hopping 
model, or XX chain. The eigenvalues follow in this case from the simple correlation 



20 
15 



-15 
-20 

















L=20 — ■ — 




I =10 




L=100 



-10 




10 



Figure 4. Size dependence of the density-matrix spectra in a critical system. Shown 
are results for segments of different lengths in an infinite hopping model. Left: Single- 
particle eigenvalues Right: Total eigenvalues w n . After Ref. [15] . 



matrix 



Cn 



" F do 

" g— iq(m— n) 



kp 



2vr 



sm(kp(m — n)) 



7r [m — n) 



(27) 



where kp = vr/2 for half filling. One sees that the whole dispersion curve is shifted 
towards the horizontal axis and becomes flatter as the length increases. The shift is not 
rapid, the first few eigenvalues vary as l/(ln L + b) with somewhat different constants b 
around 2.5. From a continuum approximation for the eigenvalue problem, one obtains 
the asymptotic formula 

r 2 



ei = ± 



7T~ 



21nL 



(21 



I 



1,2,3. 



(28) 
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which can also be derived with conformal considerations [22] • A similar expression for 
bosons was given in [53] . The formula is also valid for a segment of L sites at the end of 
a chain, if one substitutes 21nL — > ln(2L), which increases the values roughly by 2. It 
predicts the 1/lnL behaviour, but also a linear dispersion as in the non-critical case. In 
practice, this can only be seen if in addition to L also In L is large, which requires huge 
sizes. Nevertheless, it is an important guide for the understanding of the situation and 
will be used again later. Formulae of this type and the numerical difficulties in verifying 
them are known from studies of critical finite-size CTM's [HU [55l [56] . 

Although the change of the E\ is slow, it has a clear effect on the w n spectra, as seen 
on the right of the figure. The decay becomes significantly slower for larger systems, 
which means that the entanglement grows with the size. Invoking conformal results, 
one can obtain the functional form of the u^-spectrum [57l 158]. Asymptotically, ( |26|) is 
still valid, but now a ~ 1/ In L varies with the length. Therefore the DMRG method 
does not work as well in this case, although it still can handle sizes of L ~ 100. 

Finally we want to show how certain modifications of the ground state affect the 
spectra. In the previous cases, the system was always half filled, which leads to a 
symmetric spectrum (±e; appear) [29j. If the filling is varied in ( 1271) . one finds that the 
£rdispersion curve is moved up or down in a similar way as the Fermi level, see Fig. 
[5j For a completely full or empty system, which is a product state (in spin language 
all spins are up or down), the E\ are all infinite and w n becomes a Kronecker symbol, 
w„ = S n> i, as it should. 

30 I 
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Figure 5. Single-particle spectra for different ground states. Left: Variation with the 
filling. Right: Variation with the number of equal-size Fermi seas at half filling. All 
results are for a segment of L = 20 sites in an infinite hopping model. 

If the Fermi sea consists of several disconnected parts, one finds degeneracies in 
the eigenvalues, if empty and full regions in momentum space have equal size. This 
is shown in Fig. [5] on the right. It looks as if one had several independent kinds of 
particles. Effectively, the dispersion then rises only with a fraction of the slope. The 
same holds in the case of non-equal Fermi seas, where the degeneracies are washed out. 
Such a situation occurs, for example, for the ground state of the chain with an energy 
current [59J. As in the previous examples, the w n then decrease more slowly and the 
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entanglement becomes larger. 

If one modifies the hopping between the segment and the environment at one 
interface, one can interpolate continuously between a homogeneous chain and one with 
an open end. [60]. The erspectrum in this case is shown in Fig. [6] on the left. As the 
bond is weakened, a region with a steeper initial ascent appears before the curve follows 
the pattern without defect. This region can be associated with the developing free end 
and remains when the bond is cut completely. If, on the other hand, the bonds at both 
interfaces are weakened, the dispersion is shifted upwards resp. downwards as a whole 
and a gap develops. In the decoupling limit it goes to infinity and the entanglement 
vanishes. 




-5 5 10 -5 5 10 

1 1 



Figure 6. Influence of interface modifications on the single-particle spectra for 
a segment with L = 50 sites in a hopping model. Left: Modified bond at one 
end. Right: Modified bonds at both ends. The curves correspond to bond values 
t = 1; 1CT 1 ; 1(T 2 ; 1(T 3 ; 1CT 4 , from bottom to top in the right part of the fi gures. After 
Ref. [60]. 



4-2. Two dimensions 

(I) Non-critical systems 

The simplest two-dimensional system consists of a set of M uncoupled identical 
parallel chains, all divided at the same point such that the subsystem has the form of 
a half-strip [61J. This is the usual DMRG geometry. The combined RDM is then a 
product of the individual ones and TC a becomes a sum 

n a = J2^fih, (29) 

l, (J, 

where /x is the chain index. Since £/ )/L4 = Ei, the single-particle eigenvalues are simply 
M-fold degenerate. For free particles, a coupling of the chains does not change this 
situation because one can separate the system into M new independent chains by a 
Fourier transformation in the perpendicular direction [25] [62]. The index /i in (129!) then 
becomes the Fourier index q. Only e^ q will depend on q and the M-fold degenerate 
levels will become bands. 
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For coupled oscillators and an infinite half-strip, the problem can in this way be 
solved exactly by invoking the one-dimensional results. One only has to determine the 
elliptic parameter k = k(q) for each Fourier component from the coupling k x in the 



chain direction and the frequency ui 2 {q) 
u(q) 1 — k 



Uq + 2k y (l — cosg) via 



(30) 



kg. k 

Numerical results for a system of 10 chains with actually finite length are shown in 
Fig. [3 The coupling of the chains was varied and one can see nicely, how the plateaus 




Figure 7. Single-particle eigenvalues for one-half of a 10 x 10 system of coupled 
oscillators with loq = k x = 1 and different couplings k y . Reprinted with permission 
from 25j. ©2000 by the APS. 



with 10 levels develop into bands and a rather smooth, roughly linear curve results 
in the isotropic case. The initial plateaus, combined with the large freedom in the 
bosonic occupation numbers, lead to even larger plateaus in the w n -spectrum. In the 
isotropic case, one can derive an asymptotic formula as fl26l) by assuming a strictly linear 
behaviour with a slope 

: I (31) 



E t = XI 



M 



Then one finds f[2oT) with a coefficient a = X3/2tt 2 . The crucial difference is that 
A ~ 1/M depends inversely on the width, which makes the decay of the w n exceedingly 
slow for wide systems. The entanglement becomes correspondingly high. This is a 
general feature and will be taken up again in section 6. For the DMRG it means that 
the width of the strip puts a fundamental limit on its applicability. 

For subsystems in the form of L x L squares embedded in an infinite lattice, one can 
obtain similar results by solving the equation ({19]) numerically. One finds again bands 
as in Fig. [TJ but the number of states in the lowest bands is now given by (4L — 4), 
which one recognizes as the the number of boundary sites. Plotting the e% as a function 
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of the scaled index Z/(4L — 4), the results fall essentially on top of each other. This is 
the same behaviour as for the single straight boundary, where l/M enters. It is a clear 
indication that the single-particle states are associated with the interface between the 
subsystem and its surrounding, as in one dimension. 
(II) Critical systems 

In this case one finds similar features which we will exhibit for the hopping model 
on a square lattice. The isotropic half-filled model has the well-known quadratic 
Fermi surface with corners at (±7r, 0) and (0, ±7r) in momentum space. This gives 
the correlation function as the product of two one- dimensional ones as in (J2j 



C(x, y\0, 0) = 2 M<x-y)/2)MHx + y)/2) 
7r(x — y) n(x + y) 

where x and y are integers. If the model is anisotropic, the Fermi surface is more 
complicated and one momentum integration has to be done numerically. With these 
functions one can calculate the spectra for arbitrary subsystems embedded in an infinite 
lattice. For half filling, the spectra are again symmetric, i.e. the eigenvalues occur in 
pairs ±e. 

Fig. [8] shows results for L x L squares, plotted to exhibit the scaling behaviour. On 
the left, the e\ are shown as a function of the scaled index l/L. One can see low-lying 
bands which all have the same horizontal length 1 and thus contain L states. However, 
their height still varies with L. Only by plotting E\ InL they all collapse on one curve, as 
shown on the right. This demonstrates that, on the one hand, the linear size L enters as 
in the non-critical case, but that also the inverse logarithmic dependence on L found in 
one dimension remains. As a result, logarithmic corrections appear in the entanglement 
entropy, see section 6. Note also that L enters, and not 4L — 4 as before. This is most 
obvious in a band of L eigenvalues which are exactly zero (the figure shows only one-half 
of it). The latter feature is peculiar to the square and does not occur for rectangles, 
where the dispersion rises smoothly from zero. 




Figure 8. Single-particle spectra for L x L squares in an infinite planar hopping model. 
Left: £i vs. l/L. Right: e;lni vs. l/L. Only the positive eigenvalues are shown. 

The resulting spectrum of p a is shown in Fig. [9] for three relatively small systems. 
For the 4x4 square, all 2 16 eigenvalues are displayed and the s-shaped curve actually 
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reflects the symmetry of the E\ spectrum. The 4x5 system gives much smoother results 
which can be fitted well by the law fl26l) . For it, and also for the 5x5 system, the curves 
drop only to a value of about 10 -4 for n around 1000, which is to be compared with 
the one-dimensional results of Fig. HJ where this value is reached already at n ~ 100 
for L = 100. The same feature is found for other geometries [281 W7\. This shows very 
clearly the basic difference between one and two (and also higher) dimensions. 





- 5x5 - 




4x5 - 




4x4 





1 2 3 4 5 6 7 

n (xlO 4 ) 



Figure 9. Total eigenvalues w n for squares and rectangles in an infinite planar hopping 
model. For the 4x4 system, the figure shows all w n . 



5. Further aspects 

5.1. Single-particle wave functions 
(I) Chains 

The eigenfunctions associated with the Si have a particular nature. In FigJTOl they 
are shown for the smallest E[ in the case of a segment in a half-filled hopping model. 
On the left, the model is dimerized, i.e. non-critical, and one sees that the amplitude 
is concentrated near the two interfaces to the remainder and almost zero in the middle. 
This feature persists even in the homogeneous critical case seen on the right, although 
there is now a slow decay into the interior. For the highest ej, on the other hand, the 
amplitude is concentrated in the centre of the subsystem and the eigenfunction resembles 
a Gaussian. The same pattern can be seen in oscillator chains [25] [63]. [36] . It is very 
suggestive, since it means that the states which are most important in the entanglement 
are those closest to the boundary. The whole entanglement appears as a phenomenon 
taking place within a layer whose width is given by the correlation length. 

A lattice result for <pi(j) in the hopping model is only available in the case Eq = 
which occurs for odd L [52]. The wave function then is u-shaped and vanishes at every 
second site, see [61]. However, one can derive an expression in the continuum limit. 
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Figure 10. Lowest lying single-particle eigenstatcs in a dimerizcd (6 = 0.1, left) and 
a homogeneous (S — 0, right) hopping model for a segment of L = 100 sites. After 
Ref. US]. 



Putting x = j/L, it reads for a segment located between x = and x = 1 [52] 

a/x(1 — x) ^ l — a; 

Such logarithmic oscillations were found earlier in CTM studies related to the TI [56] 
and the oscillator half-chain [65]. In the latter case, which was also treated in [53], the 
square- root pref actor is absent. 

(II) Planar systems 

It is clear that the basic feature, namely the concentration near the interface, will 
also be found in two dimensions. As an example, Fig. [TT] shows the situation for 



0.012 - 




Figure 11. Squared amplitudes summed over the band of states with Ei = for a 
20 x 20 square in an infinite planar hopping model. 



the ei = states which occur for a quadratic subsystem in a planar hopping model. 
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Due to the degeneracy, the individual states are not uniquely defined and one has to 
consider all simultaneously The maxima at the boundary are clearly visible and one 
has the same u-shaped pattern as in the one-dimensional case. In addition, there is a 
slight enhancement along the diagonals. In general, the eigenf unctions have variations 
parallel to the interface which are related to the square symmetry. To bring out their 
"radial" behaviour one has to calculate the analogue of a radial distribution function. 
One then sees, apart from a small bump in the centre, a clear increase towards the 
boundary for all low-lying bands. This will be even more pronounced in a non-critical 
system. 

5.2. Nature ofH a 

The eigenfunctions presented in the previous subsection have their origin in a particular 
form of the effective Hamiltonian, which we now address. In the CTM approach it is 
possible to give an explicit expression for Ti. a . This is done by considering (|2"T|) in the 
limit of a very anisotropic system [661 [671 168] . For the TI half-chain this leads to the 
result 



where the constant C depends on A. Therefore 7i a also describes a TI chain, but an 
inhomogeneous one, with coefficients which increase linearly away from the interface. 
In the two-dimensional problem, this reflects the wedge-shaped geometry. In the 
RDM context, it suppresses the influence of sites far in the interior because 7i a enters 
exponentially into p a . This Hamiltonian can also be diagonalized directly J67J [68] and 
one recovers the result ( 1221) for ei. In the limits A — > and A — > oo, the level structure 
( |22|) can directly be read off the coefficients in ([34"]) . 

For any finite subsystem, 7i a can in principle be determined numerically. This is 
particularly simple for the homogeneous hopping model. Then the matrix elements /iy 
follow from the correlation-matrix eigenfunctions via 



The result for a segment in a chain is shown on the left of Fig. []21 The dominant 
elements are those for nearest-neighbour hopping and vary roughly parabolically. This is 
the generalization of the linear law in the semi-infinite chain to this geometry. However, 
there is also hopping to more distant neighbours, although with rapidly decreasing 
amplitude. If the segment is located at the end of a chain, one finds the same behaviour 
but with only one-half of the parabola, i.e. the hopping saturates at the free end. The 
situation for a square in a two-dimensional lattice is shown on the right of the figure. 
Going parallel to an edge, the hopping in this direction varies again parabolically. It is 
smallest close to the edge and largest halfway in between the edges. This shows that 
the inhomogeneity in 7i a always follows the same pattern. One finds it also in the XXZ 
model with A = 1/2 [69]. 




(34) 




(35) 
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Figure 12. Matrix elements in H a for a hopping model. Left: First, third and fifth 
neighbour hopping in a segment of L = 16 sites. Right: First-neighbour hopping in a 
10 x 10 square. 



In the XX chain, one can actually show that Ti a for a segment commutes with the 
operator, 

T = £ [c!c m + ct +lQ ] (36) 

i=l 

where the hopping is strictly only to the nearest neighbours and has exactly parabolic 
form [52]. Thus they have common eigenfunctions, and the result fl33|) was actually 
found from T. Also the low-lying eigenvalues are related, and it could be that in the 
limit L — > oo both operators become identical up to a factor. One cannot check that 
numerically, however, because then large E\ appear which are not accessible, see section 
3. 



5.3. Definition of a temperature 

It has been pointed out in section 3 that p a is not a true Boltzmann operator, since 
7i a differs from the Hamiltonian H, as shown above. However, if the single-particle 
excitations have the same functional form, one can bypass this argument. This is the 
case for the homogeneous hopping model [70]. Then the E\ vary linearly for large L 
according to (1281) and the same holds for the single-particle energies in H in the vicinity 
of the Fermi point. For hopping to nearest neighbours with matrix element t/2 these 
are, in the subsystem, given by 

= l WVT) < a " 1] (37) 

Therefore, one can write S\ = /3u>i with an effective temperature 

T = t7r— - (38) 

which depends on the length of the subsystem and vanishes for L — ► oo. Therefore p a can 
be regarded as a true grand canonical Boltzmann distribution for all expectation values, 
where only the small single-particle energies are important and the wave functions do 
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not play a role. This holds, for example, for the particle-number fluctuations in the 
subsystem, which vary as T L at finite temperatures. Inserting fl38l) . this is turned into 
the In L-behaviour for the segment in the chain. 



5.4- Thermal states 

Although our interest is in ground-state properties, it is instructive to see what happens 
if one calculates p a for a system at a finite temperature. This is quite easy with the 
correlation function approach, and the resulting spectra for the homogeneous half-filled 
hopping model with t = 1 are shown in Fig. [13j The steepest curve is the ground-state 
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Figure 13. Single-particle spectrum as a function of the inverse temperature for a 
segment of L = 40 sites in an infinite hopping model. 



result. As the temperature is increased it flattens, bends over and assumes the shape of 
the dispersion uj q = — cosq for the single-particle energies in H. In fact, one can write, 
expanding the Fermi function for f3 = 1/ksT <C 1 



n dq 



-iq(m—n) 



^2n 2 V H 2 ' 



which has eigenvalues in the subsystem 
and gives 

ei = -(3 cos qi 



71 



L + l 



I 1 = 1,2, 



(39) 

(40) 
(41) 
(42) 



In other words, for high temperature 
H a — > f3H Q 

which is a very plausible result. Apart from the shape of the spectrum, the essential 
point is that the level spacing is reduced from a value of order one to ~ 1/L. Such a 
situation is also found in quenches, see section 7. 
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6. Entanglement entropy 

In this section, we show how the properties of the RDM spectra seen in section 4 
translate into in the behaviour of the entanglement entropy. Due to the form of the p a 
it is given by the same expression as in statistical physics 

S = ±5>(l±e-) + £^I (43) 
l i 

where the upper (lower) sign refers to fermions(bosons). From this formula, one can 
immediately draw two general conclusions 

• The largest contributions come from small E\ (corresponding to high temperature in 
usual thermodynamics). Therefore the entropy will be particularly large in critical 
systems. For fermions its maximum value of L In 2 is reached if all £/ vanish. 

• If all ei are m-fold degenerate, the value of S is m times its value without degeneracy. 
This answers e.g. how S compares for one or two (noncritical) interfaces, or for one 
or two Fermi seas. 

(I) One dimension 

Analytical results can be given for the non-critical half-chains with the spectrum 
(T22|) . The sums then lead to elliptic integrals [60] and one obtains for fermions in the 
disordered region 



s-L 

b ~ 24 



(44) 



while for bosons the formula is 

W 4 \ _ ATi h\ T( h'\ 1 

(45) 



^"24 



k 



2 



A similar expression with an additional contribution of In 2 coming from the eigenvalue 
Eq = holds in the ordered region. Also the results for XXZ and XYZ chains [72l |4"3] 
can be brought in this form. The entropy for the anisotropic XY chain with h = 
can be written as the sum of the expressions in the ordered and the disordered region 
[4"2| [7T] . A plot of S, based on a numerical evaluation of the sums, was first shown in 
|72j . Curves for the XY model can be found in [T3| - The case of a segment in an XY 
chain was treated even before the half-chain. Using the correlation matrix and solving 
a Riemann-Hilbert problem, S was obtained as an integral over theta functions [4"8l 14*9]. 
This is equivalent to the half-chain result multiplied by two. 

In the disordered region there is little difference between fermions and bosons. The 
values of S are typically of order one or smaller, so that the corresponding ground states 
have M e fj ~ 1 — 10 states in the Schmidt decomposition. This reflects the rapid decay 
of the spectrum in Fig. [3j An exception is only the vicinity of the critical point. As 
anticipated, S becomes large there and actually diverges for this geometry. The formulae 
give for k — > 1 

S^ln^) (46) 
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where c = 1/2 for the TI model and c = 1 for the bosons. Since the correlation length 
varies as £ ~ 1/(1 — k), the logarithm is of the form ln(£/a) [72J. 
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Figure 14. Entanglement entropy for segments of different size L in a one-dimensional 
hopping model as a function of the dimerization parameter 8. The development of a 
singularity in case of vanishing dimerization is clearly visible. After Ref. |15j . 




The behaviour for a finite subsystem is shown in Fig. [TUfor segments in a dimerized 
hopping model. In this case, S no longer diverges at criticality but shows a maximum 
which becomes higher with increasing L. The size dependence at the critical point can 
be obtained in a very simple way [53J. Using the asymptotic form ( 1281) of the E\ in (J43J) 
and converting the sums into integrals gives 

2 mL r r° , , , , f°° , 

ae ln(l + exp(— e)) + / de 



S 



IT 2 



exp(e) + 1 



(47) 



L jo jo 
and since both integrals equal 7r 2 /12 one finds 

S= 1 -\nL (48) 

On the lattice, this behaviour was first found numerically [T21 [32] and then by using 
the asymptotic properties of the correlation matrices [741 175] . The general formula for 
critical chains is 

S = u-lnL + k (49) 
6 

Here c is the central charge, v the number of contact points between the (singly- 
connected) subsystem and the remainder of the chain, and k a non-universal constant 
which depends on the model parameters and the geometry. An estimate for k can be 
obtained if one replaces InL — >• InL + 2.5 in (j47p . using the scaling found for the first 
few eigenvalues. This gives k ~ 0.8 for the hopping model, whereas the correct value 
is k = 0.726. As the numerics show, the logarithmic behaviour of 5* can already be 
observed in relatively small systems, where (1281) is not yet valid. Since it holds for all 
conformally invariant models [7J [72] the formula (jUJ) is a central result. 
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The interpolation between one and two contact points via a modified bond has 
already been discussed in section 4.1. Regarding the entropy, it can be described by an 
effective central charge c e g = vc/2 in (j4"9l which varies continuously between 1/2 and 
1. The spectrum on the left of Fig. [6] then leads to the result in Fig. [151 A formula 




t 



Figure 15. Effective centrai charge for one interface defect in a hopping model as 
a function of the defect strength. The dotted curve is an analytical approximation. 
After Ref. [60]. 

for c ff based on boundary conformal theory was given in [76]. The problem was also 
generalized to the case of two coupled planes [77]. Completely inhomogeneous systems 
were studied in the form of chains with extended defects [78], gradients [79], aperiodic 
[80] and random JSTJ [82], [83] couplings. On the other hand, one can consider situations 
where the subsystem is not singly connected and thus has many contact points. For 
comb-like geometries, the leading term in S then becomes proportional to L [HI]. For 
example, if the sites of the subsystem are two lattice spacings apart, one has S = Lin 2 
in the hopping model. This is a direct consequence of ( l27j) which reduces to Cij = Sij/2 
and gives Si = for all /. Conformal results for multiple intervals are reviewed in [58] . 

(II) Two dimensions 

The influence of the interface on the spectra in two dimensions has already been 
demonstrated in section 4.2. In the entanglement entropy, it leads to the famous "area 
law" which has been the topic of many investigations, see [85] for a recent review. 
Consider, for example, the non-critical half-strip of oscillators. Each band of e^ q 
contributes an amount of order M to S which thereby becomes proportional to the 
length of the interface. Expressed differently, S is the sum of the M individual g-chain 
entropies and can be written, for large M, 

S = T «i* - M r ^ T s hq (50) 

For a square-shaped subsystem where the lowest band contains as many states as there 
are interface sites, one obtains an analogous result. 
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The argument also holds for critical systems [62j. Regarding a two-dimensional 
hopping model as a system of coupled chains, the Hamiltonian reads, for t — 1, 

"1 

2 l 



H = - 



EE 

q n 



' C n,g C n+l,c? + c i+l, g c n,g) + COS q c\ 



q°n,q 



(51) 



Thus for each q- value one has a chain with a chemical potential = cosg. This 
affects the filling but does not change the In L-behaviour of 5* , which therefore becomes 
proportional to M In L. This is still an area law, but the occurrence of the second length 
disturbs the picture somewhat. The same holds for a square L x L subsystem with the 
spectrum found in Fig. El There the number of states scales as L but the value of the 
ei as 1/lnL. Thus one finds logarithmic corrections to the area law. This was proven 
exactly by constructing bounds on 5* [861 EH EE] an d an expression for the prefactor 
was given in [87] . The problem was also investigated numerically in two [391 EH] and 
three [89] dimensions, and the presence of the logarithm traced to a finite Fermi surface 
in the system. For bosonic systems, on the other hand, no logarithmic corrections were 
found in the critical limit. 

Finally, we want to comment briefly on the largest eigenvalue w± of the RDM, which 
has a close relation to 5. It plays a role in the so-called single-copy entanglement, where 
one asks which maximally entangled state one can reach from an initial state [90] . From 
ffTTT) one sees that 

Wl = - e~ E ° (52) 



where E is the smallest eigenvalue of 7i a . This can be evaluated for the non-critical 
half-chains in the same way as 5. For example, putting S% = — lnu^, one finds in the 
bosonic case 



Si = -777 



1 

24 



In 



16k 



mi 



(53) 



In the critical limit, this diverges as S and one finds that S\ — > 5/2. The same holds for 
fermions [9T)l I9"T] . One can show that this is a general result for conformally invariant 
systems [HI M, [93] . 



7. Entanglement evolution 

In this last chapter we present results on the entanglement evolution after a change 
of the Hamiltonian Ho — > Hi. This can be treated via correlation functions as before 
and leads to interesting phenomena. The simplest case is a quench, where the change 
is instantaneous and generates a unitary time evolution \ip(t)) = e~ lHlt \ip ) . If Hi is 
also a free-particle operator, the arguments work as before |94j and the RDM has the 
exponential form (Till as in equilibrium but with a time dependent operator 

L 

n a (t)=Y,ei(t)fUt)f l (t). (54) 
i=i 
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In the case of particle conservation, the eigenvalues Ei(t) follow again from the restricted 
correlation matrix, but now taken at time t 

C id {t) = {^\c\{t) Cj {t)\^). (55) 

Therefore, one only needs to determine the time evolution of the operators Cj(t) in the 
Heisenberg picture. In the following we discuss three different situations. 



7.1. Global quench 

In a global quench, the system is modified everywhere in the same way, a situation 
which can actually be realized in optical lattices [95J. Then the initial state becomes a 
highly excited state of Hi with an extensive excess energy. 

An example which illustrates the situation very well is a hopping model which 
is initially fully dimerized (5=1) and then made homogeneous (5 = 0). The time 
evolution of the correlation matrix is then given explicitly in terms of Bessel functions 
ESI 



Cm,n(^) 2 



A J-IfA I r \ I -i£(m+ro) H m ~ n ) t /q-a 

"m,n ~r „ \ "n,m+l ' "n,m—l) t 6 ^ J m _ n ^ZtJ 



(56) 



The resulting single-particle spectra are shown on the left of Fig. [161 One sees that 
the dispersion is linear near zero and that its slope decreases with time. This leads to 
the initial increase of the entropy shown on the right of the figure. For times t 3> L/2, 
however, the Si approach a limiting curve and S saturates. The asymptotic form of the 
spectrum can be obtained from the tridiagonal correlation matrix C TO) „(oo) as in section 
5.4 

0(oc) = ~(l + cosgO, « = J^7 Z > l = l,2...L (57) 

leading to 

ei(oo) = 21ntan(®/2). (58) 

The spacing of the qi is proportional to 1/L and gives an extensive entropy S = 
L(2 In 2 — 1), a value which was also found in [94J for a similar quench in the TI model. 
An initial state where the sites are alternatingly full and empty, would even give the 
maximal possible value S = Lin 2. 

The build-up of an extensive entropy is a typical signature of global quenches. It was 
given a phenomenological description in terms of emitted pairs of quasiparticles which 
create entanglement between the subsystem and the remainder of the system [93J [58] . 
In our case these quasiparticles have maximum velocity v = 1. This simple picture also 
accounts for the "light-cone effect" [971 EE] reflected in the entropy at t ~ L/2, where 
the linear increase turns into a saturation. If one starts from an inhomogeneous state 
the increase of S can also be non-linear [79J. A closed expression for S(t) in the XY 
model was given in [99J. 

From the extensivity of S one might conjecture a relation of the quench state to 
a true thermal state. But a comparison of the spectra in Figs. fT3l and [TBI shows that, 
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Figure 16. Global quench in a hopping model, starting with a fully dimerized initial 
state. Left: Time evolution of the single-particle spectrum for a segment of L = 100 
sites. After [96]. Right: Entanglement entropy with the asymptotic value. 



apart from the linear region, they are different. A calculation of 7i a (oo) via ( |35l) shows 
that it has long-range hopping which decreases as l/\i — j\ in the interior. However, 
there are cases, where the final effective Hamiltonian resembles H. This happens e.g. if 
one starts from a chain with alternating site energies ±A. Then one finds that for large 
A the asymptotic E\ have the form ( )4T|) with (3 = 2/ A. This explains the observations 
in [lOOj . In general, the emergence of a p a {oo) after a global quench may still be viewed 
as a local thermalization and is a rather general feature of one- dimensional integrable 
systems, see e.g. [1011 1102j for a rigorous treatment. 

7.2. Local quench 

A very different behaviour is obtained if one makes sudden local changes in the system, 
for example by removing defects in a hopping model. The resulting entanglement 
evolution has been investigated for various situations and geometries [96j 11031 11041 [58] . 
We will consider here the case where a finite segment is joined to an infinite half-chain 
either on one or on both sides [104] . These two setups will be called the semi-infinite 
and the infinite geometry, respectively. 

The time evolution of the Fermi operators c n (t) is again given in terms of Bessel 
functions, and in the infinite geometry the correlation matrix reads 

C m>n (t) = i n - m Y,V~ l Jrn-At)Jn-l(t)C 3 0) . (59) 

3,1 

The double sum over all sites j, I has in this case to be evaluated numerically. In the 
semi-infinite geometry, a similar expression is obtained. 

On the left of Fig. [17] we show the low- lying single-particle spectrum for the semi- 
infinite case on a logarithmic time scale. Since the segment is initially unentangled, all 
£z(0) = oo first drop and evolve to a transient regime up to t « 2L where all but one 
relax to the stationary eigenvalues of the equilibrium chain. The remaining anomalous 
eigenvalue evolves rather slowly showing avoided crossings with the already relaxed 
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levels. The large-time behaviour is therefore characterized by a slow approach to the 
local equilibrium state. 
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Figure 17. Left: Time evolution of the lowest E;(t) for L = 40. The parameter to is 
chosen such that t — 2L gives 1 on the horizontal axis. Right: Entropy evolution in 
the rescaled plateau region for L = 60. Upper curve: Infinite geometry, lower curve: 
Semi- infinite geometry. After |104j . 



The resulting entropy evolution in the transient region is shown on the right of 
Fig. [T71 For both geometries one can see a plateau with a characteristic shape but the 
height and the length are different. The latter effect is already scaled out in the figure by 
choosing r = t/L for the infinite and r = t/2L for the semi-infinite case. Using methods 
of conformal field theory |103|. 1104] , one can derive analytical formulae for both cases 



S(t) 



v- In 
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UTlt 
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where v is the number of contact points and k u is a constant which depends on the 
geometry. These curves are indicated by the dashed lines in the figure and, apart from 
deviations at the ends of the interval, are in good agreement with the numerical data. 
For t<I, Eq. ( !60l) gives a logarithmic entropy growth in contrast to the linear increase 
in case of the global quench. If L — > oo this behaviour persists for all times. 

The emergence of the plateau region can be related to a front starting from the 
defect site and propagating with unit velocity. It becomes clearly visible if one looks at 
the eigenvectors belonging to the Ei(t) in Fig. [TTJ. The plateau ends when the front leaves 
the subsystem, which also explains the doubling of the length due to reflection in the 
semi-infinite case. In addition to these traveling fronts, which represent the maximal- 
velocity excitations, there are also more subtle signatures of the slowest ones. These are 
visible as flat parts in the evolution of the anomalous eigenvalue. 

In the above examples we have considered defects which initially cut the system into 
separate pieces. However, the behaviour is rather similar, if initially the corresponding 
bonds are only weakened. Only the height of the plateau decreases. Since S(t) is 
proportional to c in Eq. (IBTl . this decrease can be described by effective central charges 
[96] . These depend smoothly on the initial defect strength and one obtains similar curves 
as in the equilibrium situation depicted in Fig. [T5j A plateau is also found for local 
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quenches in a non-critical TI chain. The difference in this case is that it becomes flat 
and does not scale with the subsystem length [104] . In summary, for a finite subsystem 
a local quench is characterized by bursts of the entanglement: a rapid development of 
a plateau region is followed by a slow relaxation towards a local equilibrium. 

7.3. Periodic quench 

As a final example, we discuss a periodic sequence of changes H <->• H\ and its effect 
on the entanglement. The change in the Hamiltonian can be either global or local. 

In the global case, we consider again the dimerized hopping model and switch 
periodically between dimerizations ±5 [105] 1106] . This corresponds to a simple 
interchange of weak and strong bonds. The time-evolution operator up to the end 
of the n-th period reads 

[/(2m-) = U n , U = U U l = e - iH ^ e -iH^r ( 61 ) 

where r is the length of a half-period. For arbitrary times between periods n and n + 1 
one has to multiply U(2nr) by an additional unitary operator. Thus, the problem 
reduces to finding the diagonal form of U which can be done analytically by a Fourier 
transformation. It is convenient to write it as a single exponential of an average Hamilton 
operator 

U = e~^ , ff = 5>(£&-i£tfc) (62) 

q 

with Fermi operators £ 9 and n q . In the case 5—1, the single-particle energies are given 
by u g = •"fq/Zr where 

cos 7 9 = cos 2 t — sin 2 r cos q (63) 

The time evolution of the entropy is obtained again from (155|) and depicted on the 
left of Fig. [18] for several values of the dimerization 5 and fixed r = 0.47T. The overall 
behaviour is an initial, step-like increase followed by a sharp bend and a final approach 
to an asymptotic value. The steps are sharp in the fully dimerized case, but for smaller 
5 they become washed out and their height AS decreases. For general 5 and r the 
entropy displays additional slow oscillations. 

The characteristics of the time evolution can be understood from the dispersion of 
v q . For 5=1 and r = tt/2 it is strictly linear, resulting in a completely regular staircase 
with AS = 4 In 2. Thereby a segment of size L becomes maximally entangled after L/4 
periods. This case also gives an exact lattice example of the quasiparticle picture in |94j. 
In the general case, v q becomes more complicated and can have several local maxima, 
which give rise to the slow oscillations in S. 

Apart from the fine structure, the picture is similar to that of the single quench. 
Both problems become identical in the limit of very rapid switching, r — > 0. Then 
the average Hamiltonian is just the simple average H = (Hq + Hi)/2 = H and one 
recovers the quench to the homogeneous chain. However, the asymptotic entropy seems 
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Figure 18. Entropy evolution in periodic quenches for L = 40. Left: Periodically 
switched dimerization for r = OAtt and various S. After [105) . Right: Periodically 
connected chains. Upper curve: r = 5, lower curve: r = 1. 

to be always larger in the periodic case, and in general is a complicated non-monotonic 
function of r |1U5] . 

On the right of Fig. [18] we show results for a local periodic quench. Here two 
halves of an infinite hopping model are periodically connected and separated. The 
subsystem consists of the first L sites in one of the initially disconnected chains. One 
sees a characteristic difference. For a large half-period r, one has a step structure like 
in the case of the global quench, and the entropy grows linearly with the number of 
the periods. This is the result found analytically in [107] by studying a continuum 
model and taking the subsystem as one of the half-chains. For small r, however, the 
entropy curve resembles the plateau of a single local quench, with an additional fine 
structure due to the switching. In this case, the entropy grows only logarithmically. 
The interpretation is that, for slow switching, the system has enough time to recover 
and thereby the entanglement gain repeats itself after each new connection. For rapid 
switching, this is not the case, and for r — > one recovers a single quench as before. 
The transition between both regimes occurs around r = n/2. The phenomenon can also 
be seen in interacting systems [108J. 

8. Conclusion 

We have shown that the reduced density matrices of free lattice models have a 
special structure. This permits to view entanglement questions in these systems as 
thermodynamic problems and provides a very clear physical picture of the situation. 
In particular, the entanglement entropy can be understood from the character and the 
scaling behaviour of the single-particle spectra, as for conventional thermodynamical 
systems. Therefore the emphasis throughout the review was on the properties of these 
spectra. In addition to presenting them for a number of important situations, we 
also discussed the character of the corresponding eigenfunctions and of the effective 
Hamiltonian itself. Thereby the role of the interface between the two parts of the 
system entered in a natural way. From the character of the eigenfunctions in the ground 
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state problem, one can say that the entanglement "resides" mainly near the interface 
|12j . Therefore the states are rather weakly entangled in one dimension, but already 
in two dimensions this is no longer true and limits the applicability of the DMRG 
seriously. On the other hand, this role of the interface is not a general feature. Not 
only at finite temperatures, but also after global quenches, the entanglement entropy 
becomes extensive and typically the whole bulk of the subsystem is involved in the 
entanglement. On the other hand, simple local quenches only lead to logarithmic effects 
and time- dependent DMRG can be done. We have only considered quenches, but there 
are also results for continuous changes, see e.g. [1091 IH0| - which one could discuss in 
the same way as here. On the whole, time-dependent phenomena should be the area of 
further applications. Of course, the study of non-interacting systems is always combined 
with the hope that they serve as guides for more realistic ones. For the DMRG this is 
certainly the case. 
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